import numpy as np

j = np.complex(0,1)

def sQuad(a, b, c):
    if a != 0:
        disc = b*b - 4*a*c
        if np.isclose(np.abs(disc), 0, 1e-16):
            disc = 0
        if np.isclose(np.imag(disc), 0, 1e-16):
            disc = np.real(disc)
        if np.isreal(disc) and (disc >= 0):
            disc = np.sqrt(disc)
        else:
            disc = np.sqrt(np.complex(disc))
        return (-b + np.array((disc, -disc))) / (2*a)
    else:
        return np.NaN    

def sCubic(a, b, c, d):
    onethird = 1.0/3
    root3j = np.sqrt(3)*j
    u = np.array((1, (root3j-1)/2, (-root3j-1)/2))
    if a != 0:
        del0 = np.complex(b*b - 3*a*c)
        del1 = np.complex(2*b*b*b - 9*a*b*c + 27*a*a*d)
        C = np.power((del1 + np.sqrt(del1*del1-4*del0*del0*del0))/2,onethird)
        x = -(b + C*u + del0/C/u)/a/3
        x = x.real[np.isclose(x.imag,0)]
        x = x[np.abs(t-0.5) <= 0.5]
    else:
        return np.NaN